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ABSTRACT 


The  problem  of  loss  of  stability  of  marine  vehicles  under  cross  track  error  control 
in  the  presence  of  mathematical  versus  actual  system  mismatch  is  analyzed.  Emphasis  is 
placed  on  studying  the  response  of  the  system  after  initial  loss  of  stability  of  straight  line 
motion.  Center  manifold  reduction  and  integral  averaging  techniques  are  used  in  order 
to  study  the  bifurcations  to  periodic  solutions  and  stability  of  the  resulting  limit  cycles. 
Numerical  integrations  are  utilized  to  confirm  the  theoretical  results  and  to  establish 
regions  of  asymptotic  stability.  The  methods  used  in  this  work  demonstrate  the 
significance  of  nonlinear  terms  in  assessing  the  final  response  of  the  system. 
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I.  INTRODUCTION 


Accurate  path  control  of  ships  and  underwater  vehicles  along  prescribed 
geographical  paths  is  a  fundamental  problem  which  is  becoming  increasingly  important, 
particularly  as  the  missions  of  ocean  vehicles  become  more  sophisticated  with  strict 
requirements  for  performance.  In  order  for  a  control  law  to  be  able  to  perform  its 
mission  in  a  realistic  operational  scenario  it  has  to  be  robust  enough  so  that  it  can 
maintain  stability  and  accuracy  of  operations  in  the  presence  of  modeling  errors  and 
environmental  uncertainties.  The  robustness  properties  of  the  design  are  particularly 
important  due  to  the  unpredictable  nature  of  the  ocean  environment  and  the  changes  in 
the  hydrodynamic  characteristics  of  the  vehicle  during  turning,  changes  in  the  forward 
speed,  or  operations  in  proximity  to  other  objects  in  the  area.  For  these  reasons,  there 
exists  a  need  for  the  analysis  of  the  robustness  characteristics  of  a  particular  control  law 
design  and  the  establishment  of  a  rational  operational  envelope  based  on  stability  and 
performance  criteria.  Previous  studies  (Parsons  &  Cuong,  1977)  showed  that  gain 
adaptation  is  highly  desirable  due  to  changes  in  the  linearized  vehicle  hydrodynamics 
with  different  operating  conditions,  such  as  depth  under  keel.  The  resulting  adaptation 
scheme  (Parsons  &  Cuong,  1980)  required  significant  vehicle  motion  which  may  be 
undesirable  when  operating  in  restricted  water,  or  in  object  recognition  and  localization 
tasks.  Integral  control  techniques  (Parsons  &  Cuong,  1981)  proved  quite  effective,  but 
neglected  the  behavior  of  the  vehicle  which  becomes  very  important  at  low  speeds  and 
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hover  operations.  Model  based  compensators  exhibit  robust  behavior  under  conditions 
of  parameter  uncertainty  which  is  as  good  as  the  classical  linear  quadratic  regulators  for 
linear  output  feedback  systems  (Healey,  1992).  Alternatively,  sliding  mode  controllers 
exhibit  very  robust  characteristics  given  an  estimate  of  the  parameter  uncertainty  and/or 
disturbances  (Papoulias  &  Healey.  1992),  (Yoeger  &  Slotine,  1985).  Sliding  mode 
control,  however,  does  not  offer  an  infinitely  robust  design,  and  it  suffers  from  a  series 
of  bifurcation  phenomena  and  loss  of  stability  unless  proper  care  is  exercised  (Papoulias, 
1991). 

In  this  work  we  analyze  the  problem  of  loss  of  stability  of  a  path  keeping  control 
law  under  conditions  far  from  nominal.  We  assume  that  the  autopilot  has  been  designed 
based  on  a  nominal  model,  whereas  the  actual  system  is  different.  For  demonstration 
purposes  we  employ  a  linear  full  state  feedback  control  law,  but  the  methods  are  quite 
general  and  can  be  used  for  other  designs  as  well.  The  main  loss  of  stability  cases 
analyzed  here  occur  in  the  form  of  generic  bifurcations  to  periodic  solutions 
(Guckenheimer  &  Holmes,  1983).  We  use  center  manifold  reduction  techniques  and 
averaging  in  order  to  capture  the  stability  properties  of  the  resulting  limit  cycles  (Chow 
&  Mallet  Paret,  1977).  Particular  emphasis  is  placed  on  the  control  gains  as  the  primary 
bifurcation  parameters,  since  they  are  related  to  gain  margins  in  linear  control  theory 
(Friedland,  1986).  We  make  extensive  use  of  numerical  integrations  in  order  to  confirm 
the  theoretical  results.  All  computations  in  this  work  are  conducted  for  the  NPS 
autonomous  underwater  vehicle  (Bahrke.  1992).  and  all  results  are  presented  in  standard 
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dimensionless  quantities  with  respect  to  vehicle  length,  7.3  ft,  and  nominal  forward 
speed,  2  ft/sec. 
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H.  PROBLEM  FORMULATION 

A.  INTRODUCTION 

The  equations  of  motion  of  an  ocean  vehicle  in  the  horizontal  plane  are  presented 
in  this  chapter.  A  linear  feedback  control  law  is  designed  based  on  the  linearized 
equations  in  yaw,  sway,  and  rate  of  change  of  heading  angle  and  lateral  deviation  error 
to  provide  path  keeping.  Loss  of  stability  is  examined  for  small  changes  from  nominal 
in  feedback  gains,  and  system  properties. 

B.  EQUATIONS  OF  MOTION 

The  mathematical  modeling  of  a  steering  system  of  a  vehicle  consists  of  the 
nonlinear  sway  and  yaw  equations  of  motion.  Newton's  equations  of  motion  in  a  mov  ing 
coordinate  frame  fixed  at  the  ships  geometrical  center  are1 

m(v+r  +xGr) = Yr + Fvv + Yr + Fvv - 

stem  m 

-  /  C^CXv-COiv+CrldC 

bow 

l  r+mxJv+r)=N  *N  v+N  r~N  v+N.  b  +N.  6. 

Z  O'  T  V  T  V  O  S  Of, 

bow  (2) 

-  J  C^{0(y<r)\v^nd(, 

stem 


1  Equations  are  nondimensionalized  with  constant  forward  speed  u.ship  length  L.  and  the  dimensionless  time 
being  t.u/L. 
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where  only  nonzero  terms  are  kept  in  the  model  and  the  symbols  are  explained  in  the 
nomenclature.  The  cross  flow  integral  drag  terms  in  the  equations  of  motion  are  very 
important  in  low  speeds,  but  for  higher  speeds  their  effect  is  much  smaller  and  the 
steering  response  predominantly  linear  Also  for  maximum  maneuverability,  the  vehicle 
bow  rudder  is  deflected  at  the  same  amount  and  opposite  to  the  stem  rudder 
Equations  (1).  and  (2)  can  be  written  in  the  form 


(3) 


where 


r=a,.v-a,,r-2>., 6 


(4) 


(5) 


ft*,,  -  (/:-W-(*6-W  < 

Dav  =  (I.-Nr)(m-Yr)-(mxa-yr)(nr-mxG)  . 
Da2l=  (m-Yl)Nv-(mxc-Nx)Yy  , 

Da^2  =  (m-Yj(Nr-mxG)  (mxc-Nx)(Yr-m), 
Db[  -  (l.Jf}(Y6-Y6)-(mxG-YrXNi-Ni)  . 
Db2  =  (m-Y^N.-N,)  , 

D  =  (I.-Nr)(m-Yx)-{mxG-Yr)(Nr-mxc) 


The  transfer  function  between  rudder  angle  5  and  yaw  angular  velocity  r  is  obtained  from 
Equations  (3)  and  (4). 


b-yS  *  (dy  ^b  |  (2^by) 


b  S2-(au*a22)S+(analz-al2a21) 


(6) 


For  low  frequency  maneuvering  motions  this  second  order  equation  can  be 
approximated  by  expanding  in  Taylor  series  and  keeping  only  the  first  order  terms  only 
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where 


r 

6 


b 

s-a 


or 


r  =  ar~bb 


(7) 


(aua22-al2a7l)(a*2lbl  -aub2) 

K  i  -anb2)~b2{ana H  ,) 

(8) 

(^21^1  )* 

(^1 1  j^i)  ^ ;  i^2 1  ) 


Equation  (7).  Nomoto's  Equation  is  vers  useful  in  control  system  design  since  no 
sw  ay  velocity  feedback  necessary .  and  it  gives  the  fundamental  turning  performance  At 
steady  state.  Equation  (b)  turns  out  to  be 


r  =  0 

ar*bb  --  0  ^ 

,  »» 

a 

From  Equation  (9)  it  can  be  seen  that  the  relationship  between  yaw  rate  and  rudder 
control  angle  are  linear  with  a  slope  b/a  Experimental  results  show  that  the  actual 
relationship  is  not  linear  Large  increases  in  rudder  angle  tail  to  increase  the  yaw  rate 
according  to  the  amount  predicted  from  Equation  (d» 

The  linear  change  in  the  time  derivative  of  yaw  rate  due  to  rudder  angle  can  be 
augmented  by  a  nonlinear  term  of  the  form  *a,r,\  This  is  introduced  to  model  the 
appropriate  speed  loss  during  turning  The  term  ’a/  has  the  same  sign  as  'a*  The  new 
yaw  equation  is  therefore 
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(10) 


r=ar+a3r3+bb 

r 


6 

Figure  1  Difference  between  Nomoto's  equations  and  experimental  results 

The  mode!  becomes  complete  with  the  expression  of  yaw  rate  and  the  inertial 
deviation  rate  from  the  commanded  path.  Therefore  the  coupled  nonlinear  equations  of 
motion  for  the  marine  vehicle  in  the  horizontal  plane  are 

*  =  r  (11) 

r  =  ar+a/^+bb  (12) 

y  =  siniji  (13) 

This  system  of  equations  forms  the  basis  for  the  controller  design. 
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C.  CONTROL  LAW 

Equations  (11),  (12),  (13)  govern  the  steering  control  of  the  model  used  in  this 
section. The  control  law  can  be  expressed  as. 


6 


=  6M/tanh( 


where  in  the  vicinity  of  ^  =  r  =  y  =  0  we  have 


(14) 
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60  =  K^+Krr+Kyy  (15) 

6  is  the  rudder  angle  and  K*  ,Kr ,  Ky  are  the  control  gains  of  the  system.  The  linear 
control  law  is  50  .  The  rudder  angle  6  is  defined  by  a  hyperbolic  tangent  function  to 


include  the  saturation  to  our  problem  as  shown  in  Figure  (3).  Saturation  occurs  at  6Mt 
which  is  the  saturation  limit  generally  can  be  taken  as  0.4  rad.  The  linearized  form  of 
equation  of  motions  in  the  vicinity  of^  =  r  =  y=0are 


II 

(16) 

r  =ar+bb0 

(17) 

y  =  i|/ 

(18) 

4  * 

6. 


Figure  3 

These  equations  can  be  expressed  in  space  state  form  as 


X  AX 


(19) 


A  is  the  Jacobian  matrix  of  the  system  and  X  is  the  state  vector 


A  = 


0  1  0 
bK,  a+bKr  bKy 
1  0  0 


(20) 


(21) 

The  characteristic  equation  of  the  matrix  A  is 

k3-(a+bKr)k2 -bK^  k  -bKy=0  (22) 

The  desired  characteristic  equation  has  the  form 

A.3+a2A.2+a1X +a0=0  (23) 

The  gains  can  be  generated  in  terms  of  the  coefficient  of  actual  and  desired  characteristic 
equation  from  Equations  (22)  and  (23), 


K 


♦ 


(24) 


K. 


(25) 


K  = 

v 


(26) 


The  desired  characteristic  equation  can  be  written  with  respect  to  the  desired 
natural  frequencs  and  some  optimum  coefficients  The  ITAE  criterion  for  a  third  order 
equation  is  <Dort  |W2) 


U» 


S3  +  1.75co  S2 +2.15  <o|jS +  (•>*=  0 

n  n  n 


(27) 


Therefore  the  control  gains  can  be  calculated  for  a  given  natural  frequency,  as 


ocj  =  2.15a>jj 
a,  =  1.75  u). 

l  n 


(28) 


D.  LOSS  OF  STABILITY 

1.  Introduction 

The  control  law  guarantees  stability  if  all  real  parts  of  eigenvalues  of  the 
jacobian  matrix  are  negative.  Equation  (23)  is  stable  with  the  chosen  control  gains.  A 
small  perturbation  in  the  control  gains  or  coefficients  of  the  system  affects  its  stability. 
In  this  section  we  presented  the  computations  for  the  critical  value  for  change  in  the 
gams  and  the  coefficients  of  the  system. 

2.  Stability  Considerations 

For  a  system  with  a  third  order  characteristic  equation,  the  Routh  Hurwitz 
criterion  (Dorf.  1992)  requires  for  stability 

o2  a,  2  ap  (29) 

a.  Perturbation  in  Ki 

The  change  in  Kv  is  defined  by  a  coefficient  c.  where  c  is  any  real 
number  The  linearized  equations  of  the  system  for  a  small  change  in  Kv  are 


1 1 


*  =r 

f  =ar+bf>0 

y  =  ♦ 

where  the  control  law  is 

60  =  cK^+Krr+Kyy 
Therefore  the  characteristic  equation  is 

\3-(a  +  bKr)\2-cbK^X-bKy  =  0 
Equation  (32)  should  meet  the  Routh  Hurwitz  criterion  for  stability. 

( a+bKr)cbK f  *  -bKy 

The  critical  value  of  c  where  a  stability  change  occurs  in  the  system 
Equation  (34) 


Using  the  coefficients,  the  critical  c  is 

c.  =  0.2658 

cmr 


(30) 

(31) 

(32) 

(33) 

is  presented  in 

(34) 

(35) 
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Ccrilicol 


wn 

Figure  5  ccrit  versus  natural  frequency  for  Kr. 


c.  Perturbation  In  Ky 

The  results  for  changes  in  Ky  are  presented  similarly. The  characteristic 

equation,  is 


X3-(a+bKr)k2-bK*\-bcKy  =  0  (39) 

The  critical  value  of  c  is 


c  = 
Ccn/r 


ai  a2 


(40) 


625 


Ccritical 


7.525 


3.763 


0 


unstable  region 


stcbie  region 


0  0.5  1  1.5  2  2.5  3  3.5  4 

wn 

Figure  6  ccnt  versus  natural  frequency  for  Ky. 


d.  Perturbation  in  a 


After  similar  calculations  the  critical  value  of  c  is 


(a2  -ra)  a,  -  a0 


a  a, 


or 


(42) 


1.28484  co. 


(43) 
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Figure  7  ccrit  versus  natural  frequency  for  a 
e.  Perturbation  in  b 

After  similar  calculations  as  with  the  other  cases  the  characteristic 

equation  is. 


A.3+(c(a2+tf)-<7)a2~ra1A.+ra0=0 


(44) 


and  the  critical  values  are 


c 


emh 


«i 

a^cc;  +a) 


16 


(45) 


(46) 


o>n  +2. 15a 

Ccnt>  ~  2.15(1.75 a>n+a) 

Examining  Equation  (46),  we  note  that  all  the  coefficients  of  polynomial 
must  have  the  same  sign,  if  all  the  roots  are  in  the  left  hand  plane. The  necessary 
conditions  for  stability  are, 


C2:0 

c(a2  +a) -az  0 


(47) 


The  system  is  unstable  for  all  values  of  the  natural  frequency  if  c  is  less  then  zero.  When 
the  natural  frequency  is  in  the  interval  0  <  w0  <  |  a/1 .75 1  stability  changes  from  stable 
to  unstable  with  an  increase  in  c  (Figure  (8))  for  positive  values  of  c.  In  the  interval 
wn>  1 2.15a  |  stability  changes  from  unstable  to  stable  while  c  crosses  the  critical  point 
again  for  positive  values  of  c.  In  both  cases  a  pair  of  complex  conjugate  eigenvalues 
crosses  the  imaginary  axis.  Stability  is  also  lost  when  c  becomes  negative,  but  this  is 
associated  with  a  real  eigenvalue  crossing  zero.  In  this  work  we  concentrate  on  the  c  >  0 
case  since  in  applications  it  is  unlikely  that  a  change  in  the  sign  of  b  will  occur. 
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0  0-5  1  ‘.5 


Figure  8  ccrit  versus  natural  frequency  for  b  The  ctn,  goes  to  infinity  asymptotically  at 
w0  =  -a/1.75. 
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ffl.  HOPF  BIFURCATION 


A.  INTRODUCTION 

Hopf  bifurcation  is  the  simplest  bifurcation,  in  which  under  the  variation  of  a  single 
control  parameter  a  stable  focus  equilibrium  bifurcates  into  an  unstable  focus  surrounded 
by  a  growing  limit  cycle.  As  c  crosses  the  critical  value  one  pair  of  complex  conjugate 
eigenvalues  of  the  linear  system  matrix  crosses  transversely  the  imaginary  axis.  In  this 
generic  Poincare'  -  Andronov  -  Hopf  bifurcation  a  family  of  periodic  solutions  coexists 
with  the  stable  unstable  nominal  equilibrium  state.  Locally  as  e  approaches  crnt  the 
periodic  solutions  are  located  on  the  two  dimensional  Euclidean  plane  spanned  by  the 
eigenvectors  of  Jacobian  matrix  of  the  system  which  corresponds  to  the  critical  pair  of 
eigenvalues  In  this  chapter  stability 
properties  of  the  periodic  solutions  are 
established  In  order  to  establish  those 
properties  the  main  nonlinear  terms  that 
dominate  the  system  are  isolated  Center 
manifold  theory  is  used  to  reduce  the  flow 
to  a  two  dimensional  manifold  The 
methinl  of  averaging  is  applied  to  the 
reduced  system  Hopf  bifurcations  are  examined  for  small  changes  in  each  gain  anil 
svstem  dvnarnics  (coefficients 


Figure  9 


B.  THIRD  ORDER  EXPANSION  OF  THE  SYSTEM  EQUATIONS 

1.  Perturbation  in  K* 

In  this  case  the  equations  of  motions  are 

'I'  =  r 

r  =  ar+a3ri  +66  (48) 

y  =  sinty 


where 


6  = 

60  =  cK^+Krr+Kyy 

The  trivial  equilibrium  state  is  characterized  by  ^=r=y=0  Taylor  expansion  of  the 
nonlinear  terms  about  the  equilibrium  gives 


sinijf 

6 


♦  --*1*0(5) 


V 


br 


0(5) 


36 


(50) 


The  system  equations  are 


♦ 

r 

y 


=  r 

=  ar  +a?r3  +bb0 


bb' 


36; 


(51) 
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In  state  space  form  they  are  written  as 


X  =  AX+g\X) 


(52) 


where  A  is  the  linearized  system  matrix,  and  g^x)  contains  the  nonlinear  terms, 


[  o 


0 


A  = 


cbKf  a  +  bKr  bKy 


0 


0 


f  i  M  =  0 
S2M  =  "■ 


36 


sal 


3,  x  i|r 

*,00  =  ~  — 


where  the  50  3  term  is  computed  as 


^0  =  cKl^  +  Klr'+Ky 

3c2Kl(Krr  +  Kyy)+3K?r2(cK^+Kyy) 
3Kly2(cK^+Krr)+6cK,KrKy^ry 


(53) 


(54) 


(55) 


With  the  computed  gains  from  equations  (24),  (25),  (26),  Jacobian  matrix  is 


A 


'  0  1  0 
-co,  -o2  -o0 
1  0  0 


(56) 


The  eigenvalues  of  matrix  A  are  computed  at  the  bifurcation  point  ccn,  from  equation 
(34)  where  a  pair  of  complex  conjugate  roots  with  real  parts  are  obtained. 
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For  the  above 
introduced, 


X.  =J. 


a, 


dr 


X,  =  -j  — 


X3  -  c2 


system  of  equations  a  transformation  matrix 


(57) 


of  eigenvectors  can  be 


The  transformation 


T  = 


0 


0  - 


N 


-ou 


0 


a. 


2 

N  “o 


—  -l/«2 


jc  =7z 


z=T'x 


transforms  the  system  into  its  normal  coordinate  form. 


z  =  T'ATz  +  T'g\Tz) 


(58) 


(59) 


(60) 


where. 


22 


(63) 


v 


1 


0 


1 


Y 


v 


\  a- 

!  a, 

\  «„ 


(64) 


The  coordinate  z,  corresponds  to  the  eigenvalue  A,  which  is  real  and  negative 
A,  X.  denote  the  complex  conjugate  pair  with  zero  real  parts  At  the  critical  point  two 
distinct  eigenspaces  are  spanned  by  the  eigenvectors  associated  with  the  two  sets  of 
eigenvalues  These  eigenspaces  can  be  viewed  as  local  approximations  of  the  invariant 
manifolds  the  center  and  the  stable  manifold  The  system  exhibits  its  essential 
hiturcational  behavior  on  the  centre  manifold  associated  with  X,  and  A  Centre  manifold 
theory  reduces  the  flow  to  a  two  dimensional  manifold  (Hopf  bifurcation).  According  to 
centre  manifold  theory  the  coordinate  z.  is  expressed  in  terms  of  z  and  z:.  but  this 
expression  is  higher  order  Equation  (64)  can  be  written  as. 


\|t  = 


r  - 


v 


(65) 


Substituting  the  equations  (6Si  to  (60)  the  system  equations  in  normal  cixirdinates 


become 


\  a2 


'll4-!  '  124"!  -2  '  134"!  ‘ 


14 4 


(66) 


a, 


\  a2 


-•  +  r2l"l  "  r22Zl  ~2  "■  r23*-1^2 


24  4 


The  tenns  r(J  are  computed  later  from  our  equations  at  the  bifurcation  point  ccrit.  For 
\ alues  of  c  in  the  neighborhood  of  the  critical  point,  equation  (66)  becomes. 


i,  =  p  e; 


coe 


N  “2/ 


^2  ~rilZl  +r\lZ\‘~l  +  r  1 3 1  Z2  +ri4Z2 


(67) 


U)  € 


N  a2 


n,  ,3  ^  2  ^  2 

r  e^2  +  *21  M  +  *22^1  ^2  +  *23^1  ^2  +  *"24^2 


The  parameter  t  is  the  difference  between  cCfi,  and  c  or.  c=ccrit  +  e.  The  terms  /3’  and  uj' 
are  the  derivative  of  the  real  and  imaginary  parts  of  eigenvalues  with  respect  to  c 
evaluated  at  c  =  ccnt.  and  computed  from  a  perturbation  series  approach.  The  perturbation 
is  expressed  as  follows  from  equation  (34), 


c  = — 2_+e  (68) 

a,a2 

The  characteristic  equation  can  be  written  as  follows  for  a  small  change  for  Ky  in  the 
neighborhood  of  ccrit. 


+  a,X2  *(—  +  a,e)A.  +a0  =0 

a. 


(69) 
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The  solutions  of  equation  (69)  are  presented  in  equation  (70) 


Ka  =  ,  —  +cie 

\  a. 


X3  =  -a2  +a{e 


The  variables  b,  and  c,  can  be  found  by  substituting  equation  (70)  to  (69)  and  neglecting 
the  successive  powers  e2,e3. 


\  =  - 


2(a0+a2> 


a2au 


2  (a0  +  a2) 

/S'  and  w'  can  be  expressed  by  the  definition  of  the  derivative  where  e  is  zero  in  the  limit. 


b,e-0 


Equation  (72)  can  be  written  in  terms  of  the  desired  natural  frequency. 


P'  =  -0.6841  a) 
o)y  =  0.5171  o. 
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a.  Calculations  of  ru  Terms 

The  terms  r.j  are  generated  by  the  last  term  of  equation  (60),  where 
T'g3(Tz)  is  the  third  order  nonlinear  part.  From  equations  (61)  and  (54) 


T-lg\Tz)= 


0 

“f(“2  +  ao) 

«23/2(“2  +  ot0) 

0 

3/2 

a2 

0 

N  ao  N  a2 


/ 


-a. 


\a2) 


3  b  s3 

Zi  ~r^o 


1  3 


(74) 


where. 


So  = 


1  ,  (  “of  3 


63 


—  +  3 


V  a2) 


2  o 
az,Z2~3 


/  CL.  \2 


V  a2) 


2  2  3 

a  ZjZ2  +a 


k  “2  / 


:,3}  (75) 


After  expanding  the  first  row  of  equation  (74)  and  organizing  the  terms,  we  get 


a2  3 

ril  =  —  («2  +  “o) 


1 

f  “o 

3&2§L 

ia2 

(76) 


12 


'  “o  1/2  3  ,  a 

—  <“2+“0)— ~ 

Va2 


(77) 


,3  s  a 

rl3  =  (“2  +  a0> 


fr282“ 


(78) 
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r!4  =  " 


a. 


—  (a2  +  ao)(a3  + 


a- 


ar 


(79) 


Similarly  after  expanding  the  second  row  of  equation  (74)  and  organizing,  we  get 


r2i  = 


-(o2  +  a0)(- 


a 


3/2 


a- 


5/2 


3*2*L  «2 


+  j_ _ 

5/2  +  6  3/2 


(80) 


a0/  3  %  a 

r22  =  (a2*ao)-— 

“2  b 


(81) 


a 


1/2 


r23 


a 


3/2 


(a2  +  a0) 


(82) 


bz6 


sat 


r24  = 


«2  +  aQ 

a. 


3^2C 


(83) 
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b.  Averaging 

Equation  (67)  can  be  written  in  a  simpler  form  . 


=  P'ezr 


*2  = 


/ 

— 

/ 

“o 

to  e  + 

— 

\ 

“2 

— 

/ 

ao 

a)  e  + 

\ 

U2) 

Z2  +  ^,(2,,^) 


Z j  +  P  £  *^2  +  ^"2 


(84) 


where  Fi  and  F:  contain  the  third  order  expansion  terms.  The  use  of  polar  coordinates 
makes  it  possible  to  decouple  equation  (84)  .  We  use  the  polar  transformation 
z,  =  R  cos(0)  ,  z,  =  R  sin(0) 

After  some  algebra,  equation  (83)  can  be  expressed  in  terms  of  R  and  0. 

R  =  $'eR  +  Fl(R  cos(0),  Rsin(0)  )  cos(0)  +  F,  ( Rcos(Q),  R  sin(0) )  sin(0)  <85) 


RQ  = 


f 

_  \ 

ao  / 

~ 

—  +  <d  e 

a2 

(86) 


R  +  F2 (Rcos(0), Rsin(0) )  cos(0)  -  F,(  Fcos(0),  Fsin(0) ) sin(0) 

Equation  (85)  can  be  written  as  follows, 

R  =  p  'eR  +  P(6)R3 

where  P(6)  is  2 x  periodic  in  the  angular  coordinate  6, 


(87) 


P(0)  =  rucos4(0)+r12cos3(0)sin(0)+r13cos2(0)sin2(0)+r14cos(0)sin3(0) 
+r21cos3(0)sin(0)  +r22cos2(0)sin2(0)  +r2,cos(0)sm3(0)  +r24sin4(0) 


Equation  (87)  is  averaged  over  one  cycle  to  obtain  an  equation  with  constant  coefficients. 


29 


R  =  $'eR  +  KR2 


(89) 


where  K  is  defined  by. 


2n 

K  =  />(0)<f0 


(90) 


Equation  (90)  is  simplified  as  equation  (91)  after  evaluation  of  the  integral. 


K  =  ^[3rn  +  r13  +r,2  +3r24] 

O 


(91) 


Equation  (91)  can  be  expressed  in  terms  of  the  natural  frequency. 


U). 


K  =  0.452 [3 a. o)^  +  — — a,  +a?)-i^lu)  ] 
b2  6^  175  2 


(92) 


Equation  (87)  can  be  expressed  as. 


0  = 


N 


—  -u j'e-FfO)/?2 

a. 


(93) 


where. 


F(0)  =  r21cos  (0) +r22cos  (0)sin(0) +  r23cos2(0)sin‘(0) +r24cos(0)sin3(0) 
-rucos3(0)sin(0)  -rpcos2(0)sin2(0)  -r13cos(0)sin3(0)  -r,4sin4(0) 


After  averaging  equation  (93),  a  constant  coefficient  M  is  formed. 
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0  =  -^w'e+MJ?2 


where  M  is  defined  by. 


M  =  —J  F(6)d6 


or  in  terms  of  natural  frequency. 


„  n-roAa  3  r  3.9686  2.0256  1 

M  =  0.7949  co„[ - + - 


b2b2saI  v^TJ5o)^ 


(1.75ao>B-  — )(—  -a2)] 
"  1.75  1.75 


2.  Perturbation  in  Kr 

The  calculation  of  the  formulas  for  the  parameters  Kr,  Kv,  a  and  b  is  similar 
to  the  calculation  for  K+.  Therefore  in  the  following  sections  only  the  main  results  are 
presented.  The  system  equations  for  a  change  in  Kr  are 


tji  =  r 

r  =  ar+a3r3  +bb 
y  =  siniji 


6  = 


( 

,tanh  - 

i 


6o  =  K^*cKrr  +  Kyy 


The  Jacobian  matrix  is 
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A 


0  1  o' 

-a,  [(l-c)a-fa,]  -aj 
1  0  0 


(100) 


The  eigenvalues  of  the  Jacobian  matrix  are  calculated  at  the  bifurcation  point  from 
equation  (38). 

*3 

The  state  vector  is  presented  in  normal  coordinates  where  z;  corresponds  to 

♦  =  -t 


V  I 


(101) 


r  =  ~s a,  z, 

1 

V  =  — IT  ~ 

V°1 

The  system  equation  in  normal  coordinates  are. 

■-1  =  P  w  e  ♦  y  a,  )  C-,  j *riz~i  zz 
4  =  (w'e +v'a^)  P  *  rZ\z\  ~  r22::  zz 


(102) 


(103) 


where. 


a.  (a, -a) 

‘  ‘  =  < 

♦ a ) 

2(  aj  ♦  al) 


0  45428 (  1  75  w  -a) 


0  1441(1. 75 u>  -at 


(104) 


The  coefficients  of  the  third  order  expansion  terms  are. 


“oai  !  l 


cif,  *  a  1 1  3  b~  6 


rra> 


(105) 


c 2  ao  a ; 

®o*“! 


(106) 


a‘  aoai 


(107) 


a0«f?  a’ 

".4  -  —  -a, 

“o~“i  I 


(108) 


ii/:  :  ic 


«o*°i 


(109) 


b:b'  al  +  a 


(110) 


(Ill) 


9/2 


r23  = 


b2*L  “o  +  «? 


'24 


2  3 

“0  +  Ctl 


3*25Lj 


(112) 


The  cubic  coefficients  K  and  M  are  obtained  after  averaging 


^=0.052824(13.8675^^-0.5 a)n+^^-(o)^2.15aa)2)(2.15o>2^2)]  (U3) 

b2dl 


M  = 


0.01 142  [20.3337^3  w,  -0.7331  go, 
20.9727 


■(a-4.6225(i)n)(2.15o)2+a2)] 


*>2sL 


(114) 


3.  Perturbation  in  K, 

The  formulas  for  perturbations  in  Kv  are  presented  in  this  section.  The 
equations  of  motion  are. 


=  r 

i  =  ar+a3r3+6  6 
y  =  sini|r 


b  -  S^tanhj-^- 


6o  =  KtV~Krr+cKyy 


(115) 


(116) 


The  Jacobian  matrix. 
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0  1  0 

A  =  -a,  -a2  -ca0 
1  0  0 


(117) 


The  eigenvalues  at  the  bifurcation  point  , 


*1.2  =  *  J 

*2  =  a2 


(118) 


The  state  equations  in  normal  coordinates, 


=  zt 


=  22 


(119) 


The  normalized  equations  of  motion, 


z,  =  p/ezl-(a)/e+/a;)  Zj +ruz,3  ^zfz^  rl3ztzj +r14zj 
^  Zj+  P;ez2  +  r21z3  +  r22z12z2  +  r23z1z2+r24z23 


(120) 


where  /S'  and  co'  are  given  by. 


- 5 - =  0.0959  wn 

2(aj  +  a2 ) 

ana, 

- 5_J — _  =  0.1145 

2sfa'l(al  +  a22) 


The  coefficients  for  the  third  order  terms, 


(121) 


35 


r 


11 


“l°2 


2  2 

«l+«2 


l 

6 


a5  ) 


3b2b 


sat 


r 


12 


5/2 

a  gi  «2 

.  2»2  2 
b  ?>sat  al*a2 


_  a 2  «ig2 

13  -  2  2 

*  6*a  ai+a2 


r!4  = 


3/2 

«1  tt2 


a.+a. 


.3  \ 


a,  + 


3b2  b 


2*2 


a2Vai 


'21 


a2  +  <tl 


1 _ 

3i>26L  a 


6 


r22 


6^,  a,  +  a2 


5/2 


r23  ~ 


.  2  t  2  2 

*>  6*»  «l*a2 


r24 


a,  +  a. 


36z6l 


+a. 


(122) 


(123) 


(124) 


(125) 


(126) 


(127) 


(128) 


(129) 


The  cubic  coefficients  K  and  M  , 
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2  2  15  a>;  ,  ,  (130) 

K  =  0.05 1 55  [ 6.45  0.875  - -  ( 1  75  +a)(2  15  g>;  -  a*)] 

*2»L 

M  =  0.3969a,  u>(  1.75a  2.15u>  )(  2. 15  u"  *  a2) 

*  #1  *»  fi  FI  n  n 

(131) 

4.  Perturbation  in  a 

The  equations  of  motion. 


iji  =  r 

r  car  *a^rl  * bb  (132) 

y  =  sin  41 

where. 


The  Jacobian  matrix. 


6  6M,tanh|  6° 

60  =  K,*.K,r-K%y 


A 


0  1  0 

a,  a(r  I  i  a,  a(l 

l  0  0 


The  eigenvalues  at  the  bifurcation  point. 


A., 

A, 


1  }  v“i 


(133) 


(134) 


(135) 


and  the  state  equations  in  normalized  coordinates 


I 


V0!  *2 


(136) 


The  normalized  equations  of  motion  in  the  neighborhood  of  the  bifurcation  point. 


=  P  €:,-(<*)  e>v  a,)  ^  Tnz,3  rljZlzj  *rl4z? 

=  (u'e*^)  2,*  P  eZj  +^lZl3+r222l2z2+r23Zl22  +r24Z2 


(137) 


3  2 

«1  *«£ 


0.454a 


(138) 


=  -0.1441a 


«i  ~ao 


The  coefficients  of  nonlinear  terms. 


l  5  n  v 

_  /v  - 


®5‘ail  ^ 


(139) 


1  “o«i  ,  .  ,  “o  , 

-  - (,  a .  (a,  -*-a ) - ) 

b  h«u  “0*“!  V“l 


(140) 


1  aoai  1  2  ^ 

rP  — — ,  (  (a,  -a)*a,  «o  a!  a(1(a,  +a)) 

a0-a, 


(141) 


{-a,a,  - 


a0+a, 


3i26^v/^ 


((a,+a)2  af2 


-3(a2+a)2a!/2a0  +  3aja,"1/2(a2+a)-— 


(142) 


aoai  1  ai  ao 

ao+a!,  3fc26L  ao  6ai 


(143) 


,2,2  2  3 

b  6*»  go  +  gi 


a,  aft 

-4 


(144) 


,  2t2  2  3 

*  a0  +  gi 


(-(a2+a)2aj  -ao+2a,a0(a2+a)) 


(145) 


3/2  5/2 

gQgl  gl  1  gl 

V«!  fl3go  3^26L“o 


(146) 


i(o2  -a)3a^-3(a2  +  a)2aJ'za0  +  3<Xoa1ir2(a2+a)  -  — ^  >1 


The  cubic  coefficients  K  and  M  are 


K  0.03603[ -0  7331  <^+6.777903  0); 

* — 5 — (6.7779 to*  +  ( - 1 - «■  14.5725<o*)  (-u**(4.3o|  -  1 )( 1 .75  (o  +a)) 

3,2*2  1  4662u). 


0.7166 


(1.884  (o2*  1.4663  a  to,,)’)] 


(147) 


M  =0.036031-— ^  6.45  a,  o)^3 
4.3  3 

+__L_[  -21.3675c)3  +2.150)^  -7.6314co*  -  11.8787acon  +2.152a2«J)148) 

-  2. 152  co  J  ( 1.884  o)3  +  2. 151/2a  o)n) 

+ - - - (1.884co^2.151/2acoJ3]I 

4.3988  co 

n 


5.  Perturbation  in  b 

The  equations  of  motion  are, 


i|r  =  r 

r  =  ar+a3r3  +  cfe5 
y  =  sintjr 

where. 


60  =  K^iif  +Krr+Kyy 

and  the  Jacobian  is. 


The  eigenvalues  are. 


A 


0  1  0 
-ca,  ( 1  -c)a  -ca2  -caQ 
[  1  0  0 


(149) 


(150) 


(151) 
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K*  = 


“o+a“i 


(152) 


A.  =  — 

3  fi 

“l 


The  normal  coordinate  transformation  is, 


'if  =  Zj 


r  =  - 


a  +aa, 

o  i 

a2  +a 


(153) 


a2+a 
ao  +aai 


The  normalized  equations  of  motion  are, 


z,  =  o>/e  + 


«o  +a«i 
a2  +a 


^2+rUZl3  +ri2Zlh  +  ri3ZlZ2+ri4Z2 


(154) 


a0  a  CC  <  3  2  2  3 

4  =  w'«+4  -  Zl+  P'eZ2  +r21Zl+r22ZlZ2+r23ZlZ2  +r24Z2 

\  a,  +a 
'  2 


where  0'  and  <o'  are, 


P;  =  - 


ai(a0  +  aja)(a2  +a)  2.152(wn+2.15a)(1.75o)n+a) 


2(a0  +  a,a)a5  +  ao(a2+fl) 


10.995a)  +20.8767a 

n 


(155) 
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a 


0) 


2\| 


a2  +<2 
a0  +aa  j 


1- 


a0(a0+aaj)(a2  +a) 

2  2 

(o0+aa,)a1+  a0(a2+  a)) 


(156) 


=  1.075 


1.75  g>„  +  a 

n 

\J  o)n+2.1 5a 


‘  (o  +2.l5a)(1.75w  +a)\ 

j  _  ” _ _ _ 

<^(6.3725  con  + 10.9384  <z) 


The  coefficients  of  third  order  terms  are  given  in  terms  of  the  dummy  coefficient  Q  , 

a1a0(a2+a)1/2(a0+aa1)1/2 


Q  = 


aj(a2+a)  +  a?(a0+aa1) 


(157) 


rn  =  Q 


/ 

3 

' 

1 

a0+fl“i  ai 

a0  +  aa, 

a2+fl  3£26^\ 

a,  +fl 

2  // 

(158) 


„  a  3 

r„  =  Q - ax 

12  -.2*2 


bzb 


sat 


(159) 


ri3  =G 


a2a\ 


a2  +a 
ao+aai 


(160) 


t  =  £1 

* 14 


f  ao+aai'j 

(fl«,)3 

f «o+aoO 

3' 

"fl3 

\ 

a  „  +  a 

l  2  “  ; 

3^26L 

a,  +a 

\  2  ^ 

(161) 


r21  "  ^ 


«o  1  a? 

f  «o+aai'| 

6ai  3fc26^  ao 

CU  + 

l  2  M  y 

(162) 
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r22  ® 


Ccm  *1 


ao  \ 


g0+Qgi 
a2  +a 


(163) 


2  * 


r23  -  Q 


a  a 


^>26L 


(164) 


r24  =  Q 


/  a0+aa,\3/2 


lo\ 


a2  +  a 


a, 


3^26L  &l  N 


a2  +a 
“o  +  aai 


(165) 


The  cubic  coefficients  K  and  M  are  expressed  in  terms  of  the  natural  frequency  co, 


K  = 


oo. 


8(6.3725  co„+ a  co  +  2.152a) 


(a)  +2.15a)2 

[2.15a3o> — - - 1.075(co  +2. 15a) 


1.75oo.  +a 


2-153^ 


(166) 


£>26 


2 

SOI 


[2.15  oo3(o>  +2. 15a)2  +a(1.75co  +a)(l  +  2.152a2o  )J] 


M  0.2688^/(1. 75 o)n  +a)(coB  +2.15a)  ^  con+a2.15  co 


(2.9  co  +  6.6225a) 


1.75  o+a  4.3 

n 


1 


[2.15"a) 


.  3  co  +2.15a  ,  5 

4“  "  (2.15a4o>„ 


1.75  o).  +a 


-2.153acof  +2.153a3col2 


o>  +a2.15 


1.75a>  +  a 

n 


I] 
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C.  RESULTS 


Before  we  present  the  results  of  analytic  study  of  Hopf  bifurcation  we  have  to 
assess  the  effect  of  the  a,  term.  The  a3  term  plays  a  role  in  equation  (12),  and  from  the 
nature  of  the  equation  we  can  see  that  a3  has  to  have  a  negative  value.  This  is  because 
of  the  softening  spring  characteristic  of  the  steady  state  r  -  5  curve,  as  shown  in 
Figure  10. 


§ 

Figure  10  r  -5  graph  for  different  a3  values. 


In  the  previous  sections  expressions  for  /S’  and  K  were  found  for  equation  (89). 
The  value  of  K  depends  on  only  the  nonlinear  terms  and  since  the  eigenvalues  of  the 
Jacobian  matrix  cross  the  imaginary  axis  with  nonzero  speed  the  term  /S'  is  nonzero. 
Equation  (89)  has  two  steady  state  solutions,  one  at  R=0  which  corresponds  to  the 
trivial  equilibrium  solution  at  zero  and  one  at 
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(168) 


This  equilibrium  solution  corresponds  to  a  periodic  solution  (limit  cycle)  in  the  cartesian 
coordinates  z,,  z?.  From  equation  (168)  we  can  conclude  that. 

1.  If  B' >0  then, 

a.  If  K>0,  then  unstable  period  solutions  coexist  with  the  stable  equilibrium 
for  e  <  0,  and 

b.  If  K  <  0  ,  then  stable  period  solutions  coexist  with  the  unstable  equilibrium 

for  e  >  0, 

2.  If  IT  <0  then, 

a.  If  K  >  0,  then  unstable  period  solutions  coexist  with  the  stable  equilibrium 
for  e>0,  and 

b.  If  K  <  0  ,  then  stable  period  solutions  coexist  with  the  unstable  equilibrium 

for  e  <  0. 

The  stable  periodic  solutions  form  the  supercritical  Poincare '  -  Andronov  -  Hopf 
(PAH)  bifurcation,  while  the  unstable  periodic  solutions  form  the  subcritical  PAH 
bifurcation. The  period  of  the  limit  cycles  is  computed  by  substituting  equation  (168)  in 
(94) 


T= 


2tc 


ton  +  toe  +MR^ 


2n  ,,  a >'K~  P 'm  . 

—  (1 - e)  *  O(e-) 


co 


(169) 
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c 

A 


(b)  sober  it icai  bifurcation 


Figure  11  The  two  general  types  of  Hopf  bifurcation 


The  existence  and  stability  of  periodic  solutions  in  our  cases  is  examined  in  the 
following  subsections.  The  graphs  of  a*  and  0  are  presented  in  Figures  (12)  and  (13) 
We  note  that  for  b  and  for  natural  frequency  in  the  range  |a  1 .75  |  <  u>„  <  1 2. 15a  |  the 
system  does  not  have  a  complex  pair  of  eigenvalues  Hopf  bifurcation  does  not  occur  in 


4h 


that  interval 


1.  Perturbation  in 

The  graph  oi  &  for  K^,  is  shown  in  figure  (13).  13'  is  always  negative  for  K^, 
and  the  solutions  of  equation  (91)  is  presented  in  figures  (14),  (15),  (16).  K  is  less  than 
zero  for  low  natural  frequencies.  For  those  frequencies  we  get  a  supercritical  Hopf 
bifurcation  and  a  stable  periodic  solution  exists. However  when  the  natural  frequency 
increases  the  bifurcation  shifts  to  a  subcritical  Hopf  bifurcation.  Figure  (14)  is  obtained 
for  a3  =  0  and  various  5^.  Figure  (15)  is  obtained  for  a3  =  -3  and  various  osat.  Figure 
(16)  is  obtained  for  5^=0. 4  and  various  a3.When  |a3|  increases  the  domain  of  the 
natural  frequencies,  which  causes  supercritical  Hopf  bifurcation,  increases  as  well. 


Figure  14  KKi  versus  un  for  a3=  0  and  various  5sat. 
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Wn 

Figure  15  Kk^  versus  o>n  for  a3  =  -3  and  various  5^. 


1 - 1 - 1 - - - i _ i _ _ _ <L_i _ I 
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Figure  16  Kk^  versus  con  for  5iat=0.4  and  various  a3. 
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2.  Perturbation  in  K, 


The  term  0'  for  K,  changes  its  sign  at  o>0=  j  a/1 .75 1  and  can  be  observed 
from  Figure  (12).  Therefore,  for  values  of  K<0  there  exists  a  supercritical  Hopf 
bifurcation  whereas  subcritical  Hopf  bifurcation  forms  when  K  changes  its  sign.  From 
Figure  (17)  it  is  observed  that  8at  does  not  affect  the  domain  of  the  natural  frequency 
for  supercritical  Hopf  bifurcation.  However,  an  increase  in  the  ja3j  term  increases  the 
domain  of  supercritical  Hopf  bifurcation,  as  can  be  seen  from  Figures  (18)  and  (19). 


Wn 

Figure  17  KKr  versus  ojn  for  a3=0  and  various  Ssat. 


50 


Wn 

Figure  19  versus  con  for  oMl=0.4  and  various  a3. 
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3.  Perturbation  in 

For  this  case  (3 '  is  always  positive  and  for  e  >  0  there  exists  an  unstable 
equilibrium,  therefore  for  K<0  a  supercritical  Hopf  bifurcation  exists.  Similarly  for 
K  >  0  a  subcritical  Hopf  bifurcation  exists.  An  increase  in  5sat  increases  the  domain  of  the 
supercritical  Hopf  bifurcation  over  the  natural  frequency,  see  Figure  (20).  Also  an 
increase  in  j  a3 1  increases  the  domain  of  the  supercritical  Hopf  bifurcation,  see  Figures 
(21),  (22). 

-¥ 

0.2 

0.15 

0.1 

0.05 

^  0 

-0.05 

-0.1 

-0.15 

-0.2 

0  0.2  0.4  0.6  0.8  1  1.2 

Wn 

Figure  20  KKy  versus  con  for  a3=  0  and  various  5„, 
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Figure  21  KKy  versus  wn  for  a3  =  -3  and  various  5^, 


Wn 

Figure  22  KKy  versus  wn  for  5sat=0.4  and  various  a3. 
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4.  Perturbation  in  a 


0 '  is  always  positive  for  this  case.  For  t  >0  and  K  <0  there  exists  a  stable 
period  solution  coexisting  with  the  unstable  equilibrium.  When  a3=  0,  a  change  in  5ut 
does  not  affect  the  supercritical  Hopf  bifurcation  region  of  the  natural  frequency,  Figure 
(23).  If  a3  is  non  zero,  an  increase  in  SMl  increases  the  domain  of  the  natural  frequency 
which  causes  the  supercritical  Hopf  bifurcation.  An  increase  in  |a3|  also  increases  the 
domain  of  the  natural  frequencies  for  supercritical  Hopf  bifurcation,  see  Figures  (24), 
(25). 


Figure  23  Ka  versus  wD  for  a3=  0  and  various 
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5.  Perturbation  in  b 


Previously  it  was  shown  that  there  is  no  Hopf  bifurcation  in  the  inte^  val, 

1 1.75/a |  <oj„<  j 2. 15a |  K  can  be  examined  in  two  regions,  0<wn<  1 1 .75/a |  and 
<^n>  I  21 . 15a  | .  In  the  first  region  /? '  is  always  negative,  e>0,  and  K  is  always  negative. 
In  this  region  only  supercritical  Hopf  bifurcation  occurs.  Changes  in  5^  and  |  a3 1  have 
no  effect  in  this  region,  see  Figures  (26),  (27),  (28) .  For  the  second  region  /J '  <  0,  e  <  0, 
K<0  initially  and  supercritical  Hopf  bifurcations  exist.  For  higher  natural  frequencies, 
subcritical  Hopf  bifurcation  is  observed.  Figures  (29)  and  (30)  show  that  changes  in  6^, 
or  a3  have  no  significant  effect. 


wp 

Figure  26  Kb  versus  un  for  =0.4  and  various  a3for  the  first  region. 
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Figure  27  Kb  versus  cj0  for  a3=  0  and  various  for  the  first  region. 


wn 


Figure  28  Kb  versus  w0  for  a3=  -3  and  various  5U,  for  the  first  region. 
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IV.  SIMULATIONS 


To  illustrate  the  dynamics  of  the  Hopf  bifurcation  we  simulated  the  case  in  which 
there  is  a  change  in  the  gain  Ky.  We  compared  the  simulations  with  the  analytic  studies 
presented  in  previous  chapters.  In  our  case,  loss  of  stability  occurs  at  c=3.7625  from 
stable  to  unstable  while  c  increases  its  value,  as  can  be  observed  from  Figure  (6).  In  the 
following  simulations  we  use  the  typical  saturation  angle  for  the  rudder  of  0.4  radians. 
The  a3  term  is  assumed  to  be  zero. 

K  is  less  than  zero  for  the  natural  frequency  0.4  in  Figure  (20),  therefore  a 
supercritical  Hopf  bifurcation  exists  in  this  region.  The  behavior  of  this  supercritical 
Hopf  bifurcation  is  simulated  in  Figures  (31),  and  (32).  In  Figure  (31)  cK,=3.5  is 
simulated  for  two  different  initial  lateral  deviations  (  y0  =  0.05  and  y0  =  0.5). 
Regardless  of  the  initial  conditions  both  curves  converge  to  zero  for  this  value  of  c. 
There  is  a  unique  stable  steady  state  solution  for  this  value  of  c.  Figure  (32)  is  simulated 
for  K  <  0  and  oj0=  0.4  and  the  same  initial  condition  of  y0=0.5.  When  c=  3.5  the 
simulation  converges  to  zero  as  explained  in  Figure  (31).  When  c  is  increased  above  the 
critical  value  (here  c=4.0)  we  observe  that  the  solution  converges  to  a  limit  cycle.  This 
stable  limit  cycle  coexists  with  the  unstable  equilibrium  solution. 

Figures  (33)  and  (34)  are  presented  to  show  the  effects  of  a  subcritical  Hopf 
bifurcation.  For  this  purpose  w0=1.2  is  simulated  for  Figures  (33),  and  (34).  con=1.2 
is  in  the  region  where  K>0  in  Figure  (20).  The  initial  condition  for  this  case  is 
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y0~  °-05-  When  c  is  set  to  3.5  the  solution  converges  to  zero,  but  for  the  case  where 
c  is  set  to  4.0  the  simulations  converge  to  a  stable  limit  cycle.  This  is  the  case  in  Figure 
( 1 1  -b) .  The  simulation  misses  the  unstable  limit  cycles  and  converges  to  the  outer  limit 
cycle  shown  in  Figure  (11-b).  In  Figure  (34)  the  effects  of  the  initial  conditions  are 
further  illustrated.  In  Figure  (34)  for  y0=O.O5,  the  simulations  converge  to  zero,  but  for 
y0=0.5  the  simulations  converge  to  a  limit  cycle.  For  this  simulation  the  initial  condition 
was  large  enough  to  shift  the  results  to  a  limit  cycle  rather  than  the  zero  steady  state 
solution. 

From  the  simulations  we  observe  that  in  subcritical  bifurcations  the  magnitude  of 
the  stable  limit  cycle  is  greater  than  the  magnitude  of  the  limit  cycle  in  supercritical 
bifurcations.  Similar  results  can  be  obtained  for  the  other  cases  of  variation  in  parameters 
a,  b,  and  gains  and  Ky  . 
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SUPERCRITICAL  HOPF  BIFURCATION 


i  for  a)„  =  1.2  ,  y0=  0.5  for  two  simulations 
nd  c=4  (unstable  equilibrium)  . 


HOPF  BIFURCATION 


V.  SUMMARY  AND  CONCLUSIONS 


An  analytic  investigation  of  the  nonlinear  dynamic  response  characteristics  of  a 
steering  control  law  of  marine  vehicles  has  been  presented.  Bifurcation  theory  techniques 
were  utilized  in  order  to  assess  the  behavior  of  the  system  upon  initial  loss  of  stability. 
The  main  bifurcation  parameters  were  the  natural  frequency,  control  saturation  level, 
system  parameters  and  gains.  The  main  conclusions  of  this  research  can  be  summarized 
as  follows: 

1 .  There  exists  a  critical  point  for  each  gain  and  system  parameters  for  stability  of 
straight  line  motion.  The  loss  of  stability  occurs  in  the  form  of  Hopf  bifurcation  for  each 
case.  As  the  parameter  crosses  its  critical  value,  a  family  of  periodic  orbits  ("self 
oscillations")  develops. 

2.  The  desired  characteristic  equation  from  the  ITAE  criterion  has  a  faster  response 
for  higher  natural  frequency.  However,  subcritical  Hopf  bifurcation  develops  when  the 
natural  frequency  is  sufficiently  high.  If  c  is  increased  past  ccrit  the  amplitude  of  lateral 
deviation  undergoes  a  jump  and  takes  "large"  values  (  "hard  loss  of  stability  or  hard 
generation  of  limit  cycles").  For  lower  natural  frequencies,  stability  loss  occurs  as 
supercritical  where  the  limit  cycle  amplitude  grows  continuously  ("  soft  loss  of 
stability"). 

3.  The  stationary  point  where  change  of  stability  occurs  is  related  only  to  the  natural 
frequency  of  the  system  and  the  hydrodynamic  coefficient  a. 
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4.  Higher  saturation  limit  of  the  rudder  angle  increases  the  range  for  supercritical 
Hopf  bifurcation.  Theoretically,  if  a  saturation  limit  does  not  exist  all  bifurcations 
develop  as  supercritical. 

5.  The  term  a3  has  a  similar  effect  as  the  saturation  limit  of  rudder  angle.  An  increase 
in  the  |  a3 1  term  increases  the  natural  frequency  where  supercritical  Hopf  bifurcation 
exists,  but  an  increase  in  the  |a3|  term  is  not  a  desired  condition  since  vehicle  response 
to  rudder  angle  slows  down. 

6  Transitions  from  supercritical  to  subcritical  are  very  important  for  the  design  of  the 
control  law.  A  supercritical  Hopf  bifurcation  is  preferred  over  a  subcritical  Hopf 
bifurcation  since  the  subcritical  Hopf  bifurcation  may  develop  a  rapid  dynamic  jump  to 
a  new  limit  cycle. 

Some  recommendations  for  further  research  are  as  follows: 

1 .  Comparative  studies  must  be  performed  including  the  observer  dynamics.  The 
effect  of  variations  in  gains  and  system  parameters  must  be  studied  to  predict  the 
accuracy  of  the  observer. 

2  Continuation  techniques  for  periodic  solutions  must  be  performed  for  further 
stability  analysis.  The  cases  for  simultaneous  variations  of  more  than  one  terms  must  be 
studied  to  understand  the  dynamic  response  of  the  system  including  any  nonlinear  terms. 
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